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ABSTRACT 

The  application  of  the  least  squares  method,  using  C^  piecewise 
polynomials  of  order  k + m,  k > m,  q _>  m,  for  obtaining  approximations  to 
an  isolated  solution  of  a nonlinear  m-th  order  ordinary  differential  equation, 
involves  integrals  which  in  practice  need  to  be  discretized.  Using  for  this 
latter  purpose  the  k-point  Gaussian  quadrature  rule  in  each  subinterval,  the 

discrete  least  squares  schemes  obtained  are  close  to  collocation,  on  the 
same  points,  by  piecewise  polynomials  from  Cm_1. 

We  prove  here  that  under  smoothness  assumptions  similar  to  those 

made  by  de  Boor  and  Swartz  for  the  collocation  procedure,  i.e.  that  the 

solution  be  in  C , an  optimal  global  rate  of  convergence  0(  IaI^*™) 

is  obtained  in  the  uniform  norm  for  the  discrete  least  squares  schemes, 

provided  that  the  partitions  A are  quasiuniform.  In  addition,  a super- 

i 1 2k 

convergence  rate  of  0(  |A|  ) is  obtained  at  the  knots  for  those  deriva- 

tives I which  satisfy  0 < l < 2(m  - 1)  - q. 
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DISCRETE  LEAST  SQUARES  APPROXIMATIONS  FOR  ORDINARY 


DIFFERENTIAL  EQUATIONS 


Uri  Ascher 


1 . Introduction 

One  of  the  well-known  methods  for  approximately  solving  differential 
equations  using  piecewise  polynomial  functions,  is  the  least  squares  method. 
It  has  been  rigorously  treated  for  linear  elliptic  partial  differential  equations 
by  Bramble  and  Schatz  [8]  and  others,  and  "optimal"  global  rates  of 
convergence  were  established  in  the  L^-norm. 

The  application  of  the  metnod  involves  integrals,  which  in  practice 
should  almost  always  be  replaced  by  quadrature  (or  cubature)  sums.  A 
related  question  is  then,  how  this  discretization  should  be  performed  in 
order  to  retain  the  high  order  accuracy  of  the  continuous  method  without  using 
too  much  computing  work. 

In  this  paper  we  deal  with  ordinary  differential  equations. 

For  the  discretization  of  the  integrals  involved  in  the  least  squares  method 


here,  we  use  Gaussian  quadrature  rules,  an  obvious  choice.  The 
resulting  schemes  are  then  close,  in  a sense,  to  the  collocation  method 
of  de  Boor  - Swartz  [5],  In  fact,  the  functions  involved  are 
evaluated  or  exactly  the  same  Gaussian  points,  the  difference  between 


the  methods  being  in  the  degree  of  continuity  of  the  approximation  spaces: 
in  [ 5]  the  piecewise  polynomial  space  is  chosen  so  that  there  are  enough 
degrees  of  freedom  to  satisfy  the  differential  equation  exactly  on  the 
Gaussian  points,  whereas  here  the  degree  of  continuity  of  the  approximating 
functions  at  the  knots  (joints)  is  raised,  tying  up  some  free  parameters 
and  resulting  in  an  overdetermined  set  of  equations  which  is  solved  in  the 
weighted  least  squares  sense. 

We  shall  prove  here,  under  the  same  assumptions  of  smoothness  on 


I 


the  problem  as  in  [ 5],  that  an  optimal  global  rate  of  convergence  is 
achieved  in  the  uniform  norm,  when  using  the  discrete  least  squares 
procedure  outlined  above.  In  addition,  a superconvergence  result  similar 
to  that  in  [5],  but  for  less  derivatives,  will  be  proven.  These  results 
extend  those  of  [19],  which  verified  a conjecture  made  in  [18]. 

In  order  to  make  these  statements  more  precise,  we  introduce 
some  notation  at  this  point. 

Consider  the  linear  problem  of  finding  u(')  which  satisfies 

m-1 

(1. 1)  Lu(t)  = Dmu(t)  - Yj  a.(i)D  u(t)  = f(t),  a < t < b,  m>0 

f = 0 

and 

(1.  2)  Pi(u)  =0,  i = 1,  . . . , m 
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n p 


,/ 


J 


i.e.,  each  (3^  is  a linear  combination  of  evaluations  at  a and/or  b 
of  derivatives  of  order  < m - 1. 

Let  at  ( C*[  a,  b] , t = 0, . . . , m - 1,  and  assume  that  (1.1)  - (1.2) 
is  uniquely  solvable,  i.e.  that  Green's  function  G(t,  s)  for  this  prob. 
exists,  whence 


(1.4) 


where 


(1.5) 


D*u(t)  = J G^(t,  s)f(s)ds 


X ^ • • -X-  Ar 

\v  ' 


^(t,  s)  = (— ) G(t,  s),  I = 0,  . . . , m - 1 


■- 

\ t' 

\ 


We  approximate  u by  an  element  of  the  piecewise  polynomial  spline 
space  Sq(  A),  defined  as  follows.  Let  A be  a partition  of  I = [a,  b] 


/ 


r 


(1.6) 


ii;a  = t0<t1<...<tN.b 


I.  = (t.  .,  t ),  h.  = t.  - t.  ,,  i = 1,  . . . , N 

l ' 1-1’  l ’ i l 1-1 


|a|  = h = max  h.  . 
^ 1 < i < N 1 


Let  P^E)  denote  the  class  of  polynomials  of  degree  < r on  EC  I. 
Define,  for  0 < q < r - 2, 

S = S?(  A)  = {v  e Cq(I);  v|T  € P (I  ),  i = 1,  . . .,N}  . 


(1.7) 


‘I.  r'  i' 

i 


Let  now  r = k + m,  k > m,  q > m.  The  approximate  solution  will  be  sought 
in  the  space 

(1.8)  S = Sq+m(A)  = Sq+m(A)  0 {v;  (3,(v)  = 0,  l = l,...,m). 


> 
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•v  . 


The  discrete  least  squares  approximation  will  be  defined  on  Gaussian 


points:  with 


-1  < px  < • • • < pk  < 1 


the  zeros  of  Legendre  polynomial  of  degree  k,  define  the  Gaussian 


points  as 


Tij  = 2 (ti-l  + \ + PjV’  1 = • • • * N>  i = 1.  • • • , k • 


Also,  let  w^,  . . . , w^  be  the  quadrature  weights  associated  with  p^,  . . . , p^. 

Throughout,  let  K,  C be  generic  constants,  independent  of  h. 
Smoothness  assumptions:  The  following  smoothness  assumptions  will  be  in 
force  throughout:  For  some  integer  n with  0 < n < k,  we  have 


(1.10) 


f,  a,  « C(n+k)(I),  ( = 0 m - 1 


u . C(n+k+m)(I)  • 


It  was  shown  in  [ 5]  that  if  we  collocate  on  {t  } defined  in  (1.9)  by 
° m~ 1 

v c < S^+m(A),  we  obtain,  under  these  smoothness  assumptions,  a global 


accuracy  of 


(1.11) 


>'(u  - vc)IL  (I)<Chk+min(n'm-'>,  1 = 0,..., 


while  at  the  knots 


(1.12) 


! D*(u  - vc)(^)  I £ Ch*^11,  t = 0,  . . . , m - 1 . 


If  we  now  try  to  collocate  on  {t  } by  "smoother"  spaces 

Sq  (A),  q > m - 1,  we  obtain  more  conditions  (linear  equations)  than 
k+m 


-v  ;£!K9H& 


f.  ' 


t 


r 

J 

v 


1 


free  parameters  ("unknowns"),  i.e.,  an  overdetermined  set  of  linear 
equations.  We  weigh  the  equation  corresponding  to  t by  h.w. 
and  solve  the  resulting  system  via  least  squares  approximation.  This 
can  also  be  viewed  as  a discretization  of  the  integrals  resulting  from 
the  application  of  the  least  squares  method  to  solve  (1.1)-(1.2),  by  the 
Gaussian  quadrature  rule.  Precisely,  define  for  z,  v « ^(1) 


(1.13)  { 


N 

(z,  v).  = f zvdt , l<i<  N;  (z,  v)  = Z (z,  v)  = f zvdt 
1 I i=l  I 

l 


N 

= Z (z’zK  = ( z> z) 

i=l 

N 


r 


(1.14) 


K IN 

<z,v>.  = h.  y w.z(T  )v(t..),  1 < i < N;  <z,v)  = Z (z,v). 

1 1 j = i 1 1J  i = l 1 


^ ) = <z,  z)i?  1 < i < N; 


z I = <z,  z>  • 


2 i 


(Then,  for  z • v e P2k(I  ) we  have  <z,  v).  = (z,  v).,  and  if  this  holds 
for  i = 1,  . . . , N,  we  have  ( z,  v)  = (z,  v)) . For  a fixed  q,  q > m,  the 

>\t  o 

approximate  solution  provided  by  the  discrete  method,  v « S,  is  determined  by 


(1.15) 


min  jLv  - f] 
v « S 


* o 

The  solution  provided  by  the  "continuous"  method,  w e S,  is  determined  by 


(1.16) 


min  || Lw  - f | 
w € S 


# * 

We  shall  prove  the  following  results  for  both  v and  w : 
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Assume  that  the  partition  A is  quasiuniform,  i.e.  3C  such  that 

(1.17)  h/min  h£  <C 

l<i<_N 

and  C does  not  depend  on  h.  Then,  under  the  smoothness  assumptions 
(1.10),  an  optimal  global  error  estimate  is  obtained  in  the  uniform  norm, 

(1.18)  IlD'lu  - v*)|lL  (I)  <0(hk+mi"(n’m-,)),  ( = 0 m. 

Also,  at  the  knots,  with  ^ :=  2(m  - 1)  - q > 0, 

(1.19)  Id*(u  - v*)(t.)|  < 0(hk+n),  f = 0, 

If  jo.  < 0 then  no  superconvergence  will  occur  at  the  knots. 

Note  that,  with  q = m - 1,  p = m - 1 and  (1.19)  reduces  to  (1.12). 
These  superconvergence  results  are  similar  to  those  obtained  in  [11] 
for  the  H^-Galerkin  procedure  for  the  heat  equation.  See  also  [12]. 

An  outline  of  the  remainder  of  the  paper  follows.  In  Section  2,  we 
prove  that  the  numerical  schemes  (1.15)  are  well-defined  for  all  small 
enough  h.  We  then  proceed  to  estimate  distances  between  discrete 
and  continuous  inner  products.  In  Section  3,  we  treat  the  special  case 
L = Dm.  This  serves  both  for  illustration  and  for  obtaining,  later  on, 
the  desired  results  for  the  general  linear  case.  In  Section  4,  optimal 
global  results  are  proven  in  the  L^-norm  and  in  Section  5, results  (1. 18)— (1 . 19) 
are  obtained  for  the  general  linear  case. 

Section  6 treats  the  "continuous",  or  "integral",  least  squares 
approximation  (1.16),  and  the  nonlinear  case.  For  the  nonlinear  problem, 
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a Newton  process  of  linearization  and  iteration  is  applied.  The  linearized 
problem  is  solved  as  outlined  above.  Thus,  roughly  speaking,  the  method 
would  apply  (.?s  in  [5]),  provided  that  the  nonlinear  problem  is  "decently 
linearizable"  around  an  isolated  solution  and  the  starting  guess  is  close 
enough  to  that  solution.  In  Section  7,  numerical  examples  are  given, 
demonstrating  the  claims  above  and  comparing  the  discrete  least  squares 
schemes  to  the  collocation  scheme  of  de  Boor-Swartz  [5].  The  collocation 
scheme  is  found  to  be  better  (not  by  much,  though)  in  every  respect, 
except  that  the  dimension  of  the  approximation  space  gets  larger  for  colloca- 
tion, resulting  in  some  additional  storage  requirements  for  the  coefficients 
of  the  computed  solution.  This,  however,  is  not  a serious  factor  here, 
and  we  are  led  to  conclude  that  collocation  is  generally  superior  to  discrete 
least  squares  for  ordinary  differential  equations. 


2 . Discrete  and  continuous  inner  products 


Let 


(2.1) 

and 

(2.2) 


e = v - u 


R = Lv  - f = Le 


where  v t S^+^(A)  is  the  solution  of  (1.15).  By  the  Green's  function 
representation  we  have,  for  0 < / < m - 1, 

( b 

(2.  3)  D e(t)  = / Gj(t,  s)R(s)ds  = (Gi(t),R)  = [ (R,  G^t))  - (R,  G|(t)>]  +<R,  G^t)) 

<3 

Our  main  aim  in  this  section  is  to  estimate  the  fi  st  term  on  the  right 
hand  side  of  (2.  3).  Before  that  we  shall  prove  that  (1.15)  has  a unique 
solution  v (cf.  [19]). 

Let  W ( J)  and  HS(J)  denote  the  closure  of  c°°(J)  in  the  norms 


(2.4) 


and 


(2.  5) 


= £ Hdv 


Wp(J)  j = 0 


Lp(J) 


H (J)  j - 0 “2V 

respectively.  In  case  that  J = I,  it  will  be  dropped  in  the  notation. 
Define  also 


(2.6) 


H 


s,  A 


N 

y 

L 

i = l 


HS(I.) 


1/2 


-8- 


Ihe  smoothness  assumptions  (1.10)  are  sufficient  so  that 

(i)  - K such  that  when  G (t,  • ) is  considered  as  an  element  of 
C<n)[  a,  t ] x C(n)[  t,  b J , 

(2.7)  l|DJG^(t,  -)||  < K,  tel,  0<l<m-l,  j = 0 n. 

(ii)  Given  f e HP,  0 < p < k + n,  the  unique  solution  u for 
(1.1)  -(1.2),  assumed  to  exist,  lies  in  HP+m(I)  and 

(2.8)  II  Lu  ii  < C ||u  ||  , JC  I , 

HP(J)  HP+m(J) 

for  some  C depending  only  on  L.  We  also  have  that 

(2.9)  II  Lull  > C I|  u H . 

H 

Theorem  2.1. 

The  discrete  least  squares  approximation  (1.15)  is  well-defined, 
for  h small  enough. 

Proof. 

The  system  actually  solved  is  the  linear  system  of  normal  equations 

(2.10)  Aa  = b 

❖ 

for  the  coefficients  a of  the  representation  of  v in  terms  of  a specified 
basis  of  the  approximation  space  S = S^+m(A).  The  matrix  A is 
symmetric  positive  semidefinite,  and  the  solution  of  (2.10)  will  provide  the 
unique  minimum  of  (1.15)  iff  A is  positive  definite.  This  will  hold, 

o 

in  turn,  iff  for  any  v e S,  v ^ 0,  we  have  I Lv  J > 0.  We  shall 
now  prove  this  last  statement,  for  sufficiently  small  h. 
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Let  v t S.  Then  E^v  | =0  for  j > k + m,  i = 1,  . . . , N.  Thus 

i 

IJ  v II  = II  v ||  s > k + m . 

Hs,  A Hk+m-l,  A’ 

2 

Expanding  (Lv)  in  each  subinterval  I in  Taylor's  series  around,  say, 
t.,  up  to  terms  of  order  k + n,  and  using  the  Gaussian  quadrature 
precision  together  with  (2.8),  we  get 


I l|l.vK|  - |LV|^I  < Chk+n||V|i2k+n+m  ^ = Chktn||v||2 


H 


H 


k+m-1,  A 


In  each  subinterval,  v is  a polynomial.  Using  Markov's  inequality,  we  get 


(2.11) 


I II Lv  Jl ? - iLvif  I < Khn+1«vli2 

2 2 - Hm 


Now,  using  (2. 9)  and  (2. 11), 


1 2 H-  ii  2 r i - i 2 n - ii  2 , _ n ii  2 ...  n+1  n n 2 


lLv|  = ||  Lv  J|  + [ I Lv  | - IlLvIr  ] > C ||  v j 


For  h small  enough,  we  can  therefore  write 


H‘ 


m 


- Kh  vi 


H' 


m 


(2.12)  | Lv  J J > C ||  v II  2 


2 ' 


H 


m 


Q.  E.  D. 


We  now  prove  two  lemmas  about  the  residual  error,  R of  (2.  2). 
Lemma  2.1. 


(2.13) 

(2.14) 


There  exists  K such  that,  for  h > 0 small  enough, 

I R 1 2 < Khk 
llR|l2<Khk  . 
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V I * 

f 1 ■ 


l ) 


i i 

i « 


1 1 
(• 


•X  '( 


Proof. 


Let  v « S be  such  that 


II  Lv  - f II L (i)±°(hk> 


(e.g.  Galerkin  or  collocation  solution  [ 10] , [16],  [17],  [15],  or  an 


* 

approximation  to  u [2]).  Then,  since  v minimizes  (1. 15), 


I R 1 2 < I Lv  - fl2<  0(hk)  , 


obtaining  (2.13).  To  obtain  (2.14)  note  that 


llRll2<  ||  Lv*  - Lv||2  + II  Lv  - f II  2 • 


We  are  therefore  left  to  estimate  Lr|,  where 


* — ° 
r}  : = v - v « S . 


Clearly,  we  have  by  (2.13) 


M < 0(hk)  . 


By  (2.  8)  and  (2.12),  then, 


II  Lr,  ||  2 < C ||  r,  II  m < K I Lr,  | 2 < 0(h  ) . 
H 


Q.  E.  D. 


Lemma  2.2. 

There  exists  K such  that 
||Dk+jRll 


(2.15) 


'L  (I.) 

00  N i' 


< K,  i = 1,  . . .,  N,  j = o,  . . . ,n 


Proof- 


We 


follow  the  proof  of  lemma  4. 1 of  [ 5] . First  note  that,  by  (2. 14), 


(2.16)  |Die(t)|  = l(R,GJ(t))|  < HrH2  II Gi(t)  Il2  < 0(hk),  0<i<m-l 


^ ■ 


r W* 


v.  ; 
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Thus  II D e ||  /T  \ — ^(h^),  i = 0,  . . . , m -1.  Now,  since  R = Le, 

oo'  i' 

k+j  1 

D R on  I.  is  a combination  of  derivatives  of  af(-)  up  to  order  k + j 


with  derivatives  of  e up  to  order  k + j + m.  It  is  therefore  sufficient 
to  prove  the  existence  of  a constant  C such  that 

(2.17)  II  D^e  II  ,T  v < C,  i = l,  j = 0,...,m  + k + n. 

Loo'Y 


For  j < m we  have  (2.16).  For  j > m + k we  have  DJv  = 0 on  I 

i 

and  the  result  follows  from  the  smoothness  assumption  on  u.  For 

m < j < m + k - 1 let  u.  be  the  Taylor  expansion  for  u at  t = t 

i i 


up  to  terms  of  order  < m + k.  Then 

i * 

- d)£llC 

oo  i oo'  j'  ■ oo'  i 

For  the  second  term  on  the  right  hand  side  we  clearly  have 


l|0Jel[L  „)  < I(dV*  - u.)llL  (I  , + IId’(u.  - u)llL  (I 


lltffu.  - u)ll  „,<Khk+m-i, 

1 L00Ui; 

while  for  the  first  term  we  have,  by  Markov's  inequality  (note  that 

v - u.  e p.  (I.)), 
i k+m  i ’ 


ioV'-ulll  )<Kh",-1-Jl|Dra-1(v*-  u)ll  (I) 

oo  l 1 oo'  x' 


and,  since  clearly 


Dm  ^ L (I  ) - °(h1^  ’ 

oo'  y 


we  obtain 


I E^e  II  ,T  . < Ch*^+m  * \ m < j < m + k - 1 . 0-  E.  D. 

oo'  x 
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Next,  we  obtain  the  desired  estimates  for  an  interpolant  of  R. 


Lemma  2.3. 

Let  R(-)  e S *(A)  interpolate  R(')  at  the  points  {t  } of 
(1.9).  Then 

(a)  R exists 

(b)  Hr- ill  a v<o(h^),  i<  i<  n . 

(c)  (R,  z>.  = <R,  z>.,  i = 1,  . . . , N,  z e Loo(I)  . 

(d)  Assume  that 

(2.18)  IIr|Il  /n  < 0(hk)  . 

oo' 

Then 

(2.19a)  |(R,  G^t)).  - (R,  G^t)).  I < 0(hk+M),  t « I,  0 < I < m - 1 

where 

f n + 1,  t ^L 

(2.19b)  M = < 

l^min(n  + 1,  m - l),  t « I • 

(Note  that  if  t = t.  then  M = n + 1). 

m “ i 

(e)  Without  assuming  (2.18)  we  still  have  for  0 < i < m and  4j  € (I) 

satisfying  the  side  conditions  (1.  2),  that 

l(R.vp)  - (R , M I < 0(hk+m'i)  . 

Proof. 

(a),  (b),  (c)  are  trivial.  We  prove  (d):  Fix  i such  that  t 4 I.- 

Let  P t S ( 1)  approximate  G.(t,  • ) such  that 
n t 


Hp  " Gi(t)  HLoo(i.)  — °(h7)»  ^ \ ■ 

(This  is  clearly  possible,  by  (2.7  )).  Now, 

(2.20)  (R,  GJ(t)).  - < R,  G^(t)>.  = [(R,Gy(t)).  - (R,  P).]  + [ (R,  P).  - <R,P>.) 

+ [<R,P>.  - < R , G^(t)) . ] . 

For  the  second  term  on  the  right  hand  side  we  have  (R , P).  - <R  , P) . = 0, 

since  R • P«  P.  (I.)  and  n < k.  For  the  first  term  we  write 
k+n  i — 

kR.G^t)  - P).l  < Hr  It  L (I  jllc/t)  - P II L (I } < o(hk+n+1)  . 


2' V 


The  above  is  true  because 


SrII  sa^R11 

2'  l =»'  i' 

and  using  (2.18)  and  (b)  of  the  lemma  we  obtain 

II R II L (I  } < 0(hk+1/2)  . 


Similarly 

ii p - q.wHl  (i.)<°(hn+1/2)  • 

The  third  term  in  (2.20)  is  treated  exactly  as  the  first  one  above: 

|<R,P  - G^(t)>.  | < |R|^(I  }|P  - G^(t)  | ^ (I  } <0(hk+n+1)  . 

When  we  now  take  i such  that  t«  1.  we  have  G.(t,  • ) € Wm-*-1(I  ) 

l’  l ' 00  ' i' 

and  thus  m - t - 1 replaces  n in  the  argument.  The  definition  (2.19b) 
of  M thus  follows. 
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To  prove  (e)  note  that  our  result  will  not  be  local,  since  (2.18) 
is  not  assumed  to  hold.  We  proceed  as  above  without  breaking  to  sub- 

inter'Mls:  P S * f now  ipproximates  with 

ii„  M _ m - 1 . 

II P - w • Oi  n > 


For  R we  note  t 


<R,C  - P)  I 1 llRil  liv  - P li  2 


|R||  < !■’  R-  ^ + , R - R , ^ Oih^i 


and  combining  this  with  the  estimate  on  _ - P we  ^utain  tn< 
conclusion.  Similarly  for  the  term  | ( R , w - P>  | . 

We  are  now  ready  to  estimate  the  difference  between  tne  jiscicl<- 
and  continuous  inner  products  in  (2.  3). 

Theorem  2.2. 


With  our  smoothness  assumptions  and  assuming  (2.18)  to  hold, 


we  have 


(2.21)  |(R,G/t))  - <R,G^(t)>|  <0(hk+M),  1 = 0 m-1 


where 


(2.  22) 


min(n,m  - i),  t * t , for  all  i 


t = t,  for  some  i . 

l 


Proof. 


Since  R interpolates  R at  {t  } we  can  write,  fora  fixed  i, 
(2.23)  P.(S)  = R(s)  +Pi(s)[T.1,  ...,xik,s)R,  S e ^ 


bL  ^ 


V 


with  [ t '.  . . , t s ] R the  k-th  divided  difference  of  R and 

k 

(2-24)  P (s)  = 7T(S  - T ) . 

j=l  J 

Now, 

<R,  - <R,G,(t)>l  = (R,G,(t))(-<R,G<(t)>1  = [(i,G/(t))l  - <R,  + 

+ <P1(’)1T11 Tlk--|R.  ■ 

The  first  term  on  the  right  was  estimated  in  Lemma  2.  3(d).  For  the 
second  term  we  proceed  exactly  as  in  [ 5,  proof  of  Thm.  4.1,  pp.  600-601], 
using  Lemma  2.  2.  Now,  to  obtain  (2. 21)  we  sum  up  on  i,  i = 1,  . . . , N. 

Since  we  have  N - 1 "regular"  intervals  with  M = n + 1 and  one 
"special"  interval  with  M = min(n  + 1,  m - t),  we  obtain  an  overall 
error  (2.21)  with  M given  by  (2. 22)  for  t * t.  all  i. 

When  t = t.  for  some  i,  we  simply  do  not  have  a "special" 
interval  I.  such  that  t t I . Thus,  summing  up  N terms  of  order 
0(h*C+n+S  gives  the  order  0(h^+n).  Q.E.D. 

Corollary  2.1. 

m — / 

With  4j(-)  e W^  (I)  replacing  G^(t,  •)  in  Theorem  2.2  we 

obtain 

(2.25)  l(R,*)  - <R,  4»>  I <0(hk+m"J),  0 < / < m . 

(Note  that  (2.18)  is  not  assumed  to  hold. ) 

Proof. 

(R>  4>)  ~ (R,  4^)  = [(R»4>)  _ (R,+  )J  + [ (R  ,4* ) * ^ R > ] . 
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The  second  term  on  the  right  was  satisfactorily  estimated  in  Lemma  2.  3(e). 
For  the  first  term  we  again  write 

( R - R ,4j).  = (P(-)X')(t.1,.  • T.k,  ’ ] R) 

P.  :i/en  by  (2.24),  and  proceed  exactly  as  i [ 5,  pp.  600-601],  Note 
that  the  derivatives  of  ^ and  R are  taken  to  be  bounded  in  the  uniform 
norm,  according  to  the  assumption  on  and  Lemma  2.2.  Q.  E 


3.  The  case  of  L = Dm 

In  order  to  be  able  to  use  Theorem  2.  2 and  to  proceed  further,  we 
have  to  bound  R appropriately  in  the  uniform  norm  (i.e.,  show  that 
(2.18)  holds).  For  this  end  and  for  illustrative  purposes  we  obtain  all 
the  results  first  for  the  special  operator  L = Dm. 

Let  u be  the  solution  of  (1. 1)-(1.  2) . Then 

m-1 

(3.1)  Dmu(-)  = f(.)  + £ a (-)D  u(-)  =:  f(  • ) . 

1 = 0 

* 1c  +n 

We  have  f « C (I).  Let  U(>)  be  the  solution  of 

(3.2)  min  I Dmv  - f I ^ • 

Since  we  have  assumed  (1.  3),  we  are  able  to  conclude  that  Green's 
function  H(t,  s)  for  Dm  exists: 

/ rb 

(3.  3)  D u(t)  = J H^(t,  s)f(s)ds,  0 < l < m - 1 

a 

where 

(3.4)  H^(t,  s)  = l~J  H(t,  s)  . 

Let 

(3.  5)  v = U - u . 

First  we  prove  the  equivalent  of  (2.18)  for  this  special  case. 

Lemma  3.1. 

(3.6)  ||Dmv  ||  m < 0(hk)  . 
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V. 


\ ^ 

■ 


Proof. 


Note  that 


m ,^mTr  c 
D v = D U - f 


i.e.,  Dmv  is  the  quantity  minimized  in  (3.2).  Let  now  P « (A) 


interpolate  f at  {t.. },  i = 1,  . . . , N,  j = 1,  . . . ,k.  Then,  for  any  v e Sk+m(A), 


(3.7) 


|Dmv-fl2=  lDmv-Pl2=  II  Dmv  - P JJ  2 


(3.8) 


lf'P|lL00(I)£O(h  ) 


By  (3. 7)  and  (3.  2)  we  have  U determined  by 


D v - P 


v‘ W4) 


But  then  DU  is  j 


mU  is  just  the  L - projection  of  P into  m(A).  By  [10],  [4] 


we  have 


DmU-  P II L (I)  < 0(h  ) . 


Combining  this  with  (3.8)  we  obtain  (3.6). 


Q.  E.  D. 


We  are  now  ready  for  the  general  theorem  for  this  special  operator. 


Theorem  3.1. 


Under  our  smoothness  assumptions  (and  for  quasiuniform  partitions)  we 


have 


(3.9) 


I <o(hk+mIn(n'm-')),  1 = 0,.... 


and,  at  the  knots, 


(3.10) 


I D/v(t.)  I < 0(hk+n),  1=0, , M. 

where  ^ = 2(m  - 1)  - q.  (Note  that  for  q sufficiently  large,  < 0 and 

(3.10)  occurs  for  no  derivative,  i.e.,  no  superconvergence). 

Proof. 

For  the  pointwise  error,  we  have 

(3.11)  Div(t)  = (Dmv(-),H/(t,  •))  = [(Dmv,H/t))  - (D^.H^t))]  + <Dmv,  H/t))  . 

For  the  first  term  on  the  right,  Theorem  2.  2 applies,  since  we  have  proved 
(3.6).  For  the  remaining  discrete  term  in  (3.11)  we  have,  since  U 
solves  (3.2), 

(3.12)  |(Dmv,Hi(t))l  = l<Dmv,  H^(t)- Dmv>|  < ||Dmv  HL  | H/ (t)  - DmV  | { 

for  all  v t (A) 
k+m 

where  for  a function  ^ « 1.^(1) 

N k 

(3.13)  Ul,  :=  Z hi  L w U(r  )|  . 

1 i=l  j=l  1 J 

Now,  with  this  simple  operator,  Dmv  is  just  an  approximant  out  of  m(A) 

to  H.(t,  • ) and  can  clearly  be  chosen  so  that  [ 2] 

£ 

||Dmv  - Hf(t)||L  (I ^ < 0(hM) 


where,  if  t * I.  , 


M = 


min(n,  m - / - 1),  i^  - k < i < i^  + k 
n otherwise  . 
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Multiplying  each  of  these  estimates  by  h.  and  summing  up  on  i,  we 
obtain 


(3.14) 


mf  lH  (t)-Dmvl  <o(hmln(n'm-'))  . 

O A 


v ‘ Su  J.  ( 

k+m 


Substituting  this  and  ( 3.  6)  into  (3.12),  we  obtain  (3.9)  (via  (3.11)),  for 


0 < t < m - 1.  The  conclusion  for  l - m now  follows  from  (3.  6). 
If  t = t for  some  j,  1 < j < N,  then  we  simply  have  here 

h/<v  * sr''2(i>  • 

If  m-4-2<q-m  then  the  above  argument  can  be  applied  with 


I.  U replacing  I.  , and  (3.9)  follows.  Further,  if  our  approximation 

space  is  not  too  smooth,  then  we  can  obtain  more:  if  0 < f < ^ then 


( < 2(m  - 1)  - q implies  m-i  - 2>q-m 


and  now  we  simply  have,  instead  of  (3.14), 


min  | H .(t)  - D v|,  = 0. 
l 1 

°Q 
V € S, 

k+m 


The  result  (3.10)  follows  from  a similar  result  in  Theorem  2.2. 


4.  An  error  bound  in  the  L -norm 


In  order  to  obtain  pointwise  error  estimates  we  first  obtain  optimal 


global  results  in  the  L -norm 


Theorem  4 . 1 


k+min(n,  m-i ) 


Proof 


We  apply  a duality  argument.  Let  4*  satisfy  (1.  2)  and  solve 


where  g « W 


1,  0 < 1 < m . Assume,  without  loss  of 


generality,  that  n > m - f . Then 


(g,e)  = (L  i\j,e)  = f L 4*(t)  f G(t,  s)R(s)dsdt  = J R(s)  f G(t,  s)L  i}j(t)dtds 


We  now  make  use  of  Corollary  2.1 


For  the  remaining  discrete  term  we  have  by  (1.15) 


Taking  v to  be,  e.g.,  the  Galerkin  or  collocation  solution  to 

Ltjj  - ijj 

4>  satisfies  (1.2)  , 


we  have 


and 


Thus, 


Corollary  4.1. 
(4.4) 


!Lv  - 4J  < 0(hm_i)||Dm"i4, 


L (I) 

oc v ' 


J<R,  +)  I < IrI.Ilv  - 4,1  < o(hk+m_i)  . 


t = sup  |(g,  e)  I < 0(hKfm-i)  . 


WA(I)  g 


IrJ  ||  k+min(n,  m-i ), 

1°  e|lLi(l)  - °(h  ).  1 = 0,  • ■•-m  . 


Q.  E.  D. 


Remark. 

Note  that  the  quasiuniformity  assumption  (1.17)  was  not  used  here. 
It  is  used  in  obtaining  a bound  similar  to  (4.4)  in  the  L -norm. 
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5.  estimates  and  superconvergence  results  for  a general  linear  operator 


Here  we  extend  the  results  obtained  above  for  v,  to  e.  First 
we  establish  (2.18).  Let 
(5.1)  t|  = v'"  - U . 

Then  e = v + g and  we  are  left  to  estimate  g.  Note  that  g t S^+^(A). 
Lemma  5.1. 


(a) 

(5.2) 


Id', II,  ,I.<o(hk+mln(n’ra-')),  o<,  <m  . 


■(b)  If 

(5.  3) 

then 

(5.4) 


, 1 2 2k+l 

I Lr|  1 2 < 0(h  ) 


IrIIl  (I)i°(hki  ■ 

oc ' ' 


Proof. 


(a)  follows  immediately  from  (3.9)  and  (4.4).  To  prove  (b)  assume 
(5.  3)  to  hold.  Similarly  to  (2.16),  we  have 


(5.5) 


D t|  ||  L ^ < 0(h  ),  l = 0,  . . . , m - 1 . 


Now,  we  have  assumed 


(5.6) 


|Lt](t  ) | <0(hk),  i = 1,  . . . ,N,  j = 1,  ....  k . 


But,  by  (5.5)  this  means 
(5.7) 


DmT!(T  )|  < Q(hk),  i = 1,  . . . , N,  j = 1,  . . . , k . 


U 


Fix  i.  On  I.,  D q « p^(I.)  and  the  equivalence  of  norms  on  a finite 
dimensional  space  yields 

llDr%JI  , , < 0(hk),  i = l N. 

00 1 i* 


Combining  this  with  (5.  5),  we  obtain 


II  Lt|  ||  L 


(l)<°'h 


and  (5.4)  follows,  since  R = Le  = Lv  Lr|. 
Lemma  5.  2. 


0.  E.  D. 


(5.8) 
Prop  f . 


Under  our  assumptions  (of  smoothness  and  quasiuniform  mesh) 

Hr II,  m < 0(hk)  . 

oc''  ' 


By  the  previous  lemma  we  are  left  to  prove  (5.  3).  Since  q € S^+m(A) 


we  have 


(R,  Lq)  = 0 


i.  e. 


( Lv  + Lq,  Lq)  = | Lq  I + (Lv , Lq)  • 


m-1 


(5.9)  | Lq  | = - (Lv,  Lq)  = -(Dmv  , Dmq)  + ( Dmq,  Jj  a D V ) 


m-1 


m-1 


1 = 0 
m-1 


+ ( 


D%,Dmv>  - ( 2,  a.Dfv,  2,  a D%>  . 


t =0 


e = o 


( = 0 
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Now,  by  definition  of  v,(D  v,D  r^)  = 0.  The  other  terms  on  the  right 


hand  side  of  (5.9)  are  of  lower  order.  We  would  like  to  show  that  each 


of  the  discrete  inner  products  appearing  in  (5.9)  is  < 0(h 


is  immediate  for  all  terms  except  for  (D  r),  a 


We  show  it  for  the  latter  term 


since  on  each  subinterval  I.,  D 


estimate 


For  the  term  (D  ri,  a ,D 

17  m _ 1 


the  estimate  is  similar 


We  have  shown  (5.  3)  and  the  conclusion  (5.8)  follows 
We  are  now  ready  for  the  general  results. 


Theorem  5.1. 

Under  our  assumptions,  we  have  (1 . 18 ) and  (1.19) 

(5.14) 

and  for  0 < i < N 


D*e 


L (I) 
00  ' ' 


£0(hk+mtn(n,n >-/),  1 = 0 m 


(5.15)  |D'e(t,)l  < 0(hk+n),  1 = 0 2(m  - 1)  - q . 


Proof. 

Let  0 < i < m - 1.  Assume,  without  loss  of  generality,  that  n > m - t . 

o 

We  can  write  for  any  v e S 

(5.16)  D*e(t)  = (R.G^t))  = [(R.G^t))  - < R,  G^(t)>  ] + 

+ [<R,  G^t)  + Lv>  - (R,G|(t)  + Lv)]  + (R,  G^(t)  + Lv)  . 

The  first  two  terms  on  the  right  hand  side  of  (5.16)  can  be  estimated 
by  Theorem  2.2,  since  we  have  proven  (2.18).  We  are  left  with 


(R,  Lv  + G^(t))  . 

Now, 

I (R,  Lv  + Gi(t))  | < IIrII  JlG^t)  + Lv  || x , 
so  that  we  have  left  to  show  that 

(5.17)  mino  IlG^t)  + Lv  ||  < 0(hm"  V 

v « S 

Let  t « I . We  can  replace  G (t,  ■)  in  I by  a polynomial  piece, 
0 l0 
say,  such  that  the  resulting  function  4j  will  satisfy 

icyo  - ,nL  (I  }<  o(hm"i-1) , 
xo 


-27- 


K 

. i 


so  that 


N 


(5.18)  110/0-4.11^)=  llGi(0-4.DLi(I>)<hiollGi(t)-4.llL  (I  }<0(h 


m-i , 


and 


[D 


m-t 


0(1) 


i * i. 


L«(Ii)  lo(h_1),  i = iQ  • 


Now,  for  ijj,  clearly  there  exists  a v e S such  that 

m-i,  ii  m-t  n m-t 


(5.19)  II Lv  - +|IL  (I)  < 0(h  ) ||  D " 4i|lL  (I)  < 0(h 


) • 


Collecting  ( 5. 18)  and  ( 5. 19)  we  obtain  ( 5. 17),  hence  (5.14),  for  0<l<m-l. 
The  result  for  / = m now  follows  from  (5.8). 


As  in  Theorem  3.1, we  obtain  nothing  better  for  t = t and 


t > 2(m  - 1)  - q.  We  now  prove  the  superconvergence  result,  assuming 
n > m - t , with  t < 2(m  - 1)  - q. 


Let  now  t = t.  for  some  i.  We  break  the  interval  [ a,  b] 


into  two. 

In  [a,t  ] we  have  G (t.,  • ) « W^[  a,  t.  ] and  in  [t.,b]  we 
1 * 1 1 1 


have  G (t  , • ) « W^[  t , b] . Also  G/t  , • ) e Cm  * 2[a,b]  and,  since 


t < 2(m  - 1)  - q,  we  get  m-i-2>q-m,  and 


G/t.,  ■)  € Cq~m[  a,  b]  . 


Let  now  <)>(•)  satisfy  (1.  2)  and 
(5.20) 

Such  <j>  exists  by  our  basic  assumption  on  L and,  moreover,  we  have 


L<)>( ’ ) = G/t.,  • ) on  I 
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But  then, 
of  4>  at 

(cf.  [ 10)) 

(see  e. g. 
and  (5. 15) 


W 

W 


n 

00 

n 

00 


<J>  « Cq[  a,  b]  . 

we  can  find  Vj  € Sq+m(A)  to  "fill  the  hole"  of  low  continuity 

o 

i.e.,  3v^  e S such  that 

v2  + <j>  £ Wn+m[  a,  b] 


- Now,  for  v1  + 4>,  clearly  there  exists  v £ S such  that 
II  Dj(v  - (V]  + $))  ||L  ^ < 0(hn+m-j),  0 < j < m 

[2]).  Denoting  v = v - v we  obtain 

llLvtG/(V*l  (1)<°(hn> 

00 v ' 

foilows . Q.E.D. 
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6.  Continuous  least  squares  and  nonlinear  problems 

In  this  section  we  extend  our  results  to  the  integral  least  squares 
method  and  to  nonlinear  problems;  first,  the  integral  (or  "continuous") 


least  squares  method. 

Our  effort  was  concentrated,  up  to  here,  on  the  discrete  formulation 
(1.15),  because  this  is  the  approximation  scheme  actually  calculated  with. 
Note,  though,  that  by  a similar,  and  shorter,  analysis  the  results  are 
immediately  extendible  to  the  continuous  least  squares  method  (1.16). 
Theorem  6.1. 

jjc 

Let  w be  the  solution  of  (1.16).  Under  our  smoothness  and 

quasiuniformity  assumptions,  results  similar  to  (1.18),  (1.19)  are  obtained 

* * 

for  w replacing  v . 

Proof. 

We  simply  go  through  the  much  simplified  versions  of  the  proofs 
for  the  discrete  case  in  Sections  3,  4,  5,  and  apply  them  to  the  continuous 

case.  Q,E 

We  now  consider  the  extension  of  our  results  to  the  nonlinear  ordinary 

differential  equation 

(6.1)  . Dmu(t)  = F(t,  u(t),Du(t),  ...,Dm_1u(t)),  a<t<b 

subject  to  the  side  conditions  (1.  2).  (If  the  side  conditions  were  non- 
homogeneous,  a reduction  to  the  homogeneous  case  is  done  in  a trivial 
manner).  We  proceed  in  the  same  way  as  in  [ 5],  [17],  i.e.  apply  a 
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local  linearization  and  iteration  process  (Newton's  method).  Since  the 
method  of  linearization  is  independent  of  the  numerical  method  used  to 


approximate  the  linearized  problem,  the  results  of  [ 5],  [17]  will  apply  here. 

To  be  able  to  use  any  of  these  numerical  methods,  the  sought 

solution  u of  ( 6 . 1)— (1 . 2)  needs  to  be  an  isolated  one.  This  will  happen 

if  the  linearized  problem  at  u is  uniquely  solvable:  Consider  the 
m +1 

curve  C C ]R  defined  by 

C = {[  t,  u(t) , . . . , Dm  Xu(t) ]T,  t « [ a,  b]  } 

and  suppose  that  F t C*[7?  ],  where  V C IRm+*  is  some  6-neighborhood 
of  C.  Define,  for  0 < f < m - 1, 


(6.2)  Fi(t;v)  (t,z0 zm_1)| 

* z.  = D v,  i = 0,  . . . ,m-l 


, v t Cm  1 [a,  b] 


(6.3) 


m-1  . 

L = °m  - t a.(->D 


(6.4) 


ai(t)  :=  Fi(t;u),  i = 0,  . . . , m - 1 


Then,  if  Green's  function  G(t,  s)  for  L of  (6.3)  exists,  implying  unique 
existence  for  the  linearized  problem,  this  is  sufficient  to  guarantee  that 
3(7  > 0 such  that  u is  the  unique  solution  of  (6.1)-(1.  2)  in  the  sphere 
B(Dmu,<r)  [17]. 


For  Newton's  method  to  converge  quadratically  we  have  to  assume  more 
about  F.  A sufficient  assumption  is  that  F t C ^(v).  The  actual 
process  is  as  follows:  a sequence  of  approximations  ^vr^r>0  is 

o ° 

determined  with  « S being  a (good)  initial  guess  and  v , « S being 
o r+i 


the  discrete  least  squares  solution  to 

(6.  5)  min 

v « (A) 

k+m 


Lv-f| 
r r 2 


where 


m-1 

L v(t)  = Dmv(t)  - Yj  F ( t;  v )D  v(t) 
r 1 = 0 


(6.6) 


f(t)  = F(t,  v . . . , Dm_1v  ) - £ F (t;v  )D%  (t) 

r r * • _ r\  * * 1 


m-1 


1=0 


„m 


For  purposes  of  analysis  it  is  more  convenient  to  consider  y = D u, 
for  which  (6. 1)-(1.  2)  are  written  as 
(6.7)  y = Ty 

where 


(6.8) 


(Tf)(t)  = F(t,  f(”m)(t) f(_1)(t)) 

f(_m+i)(t)  = (H^(t),  f),  / = 0,  . . .,m  - 1 , 


H (t,  s)  given  by  ( 3.  4).  The  approximation  Y s S/1  m(A)  (to  which 


(D  v } >Q  will  eventually  converge)  that  we  wish  to  find  satisfies 


(6.9) 


Y = P^TY 


where  is  the  discrete  L^-projector;  i.  e.  P^  associates  with  a 
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function  ft  C[a,b]  an  approximation  P^f  i m(Y)  such  that 


(6.10) 


|PAf  - f 1 2 < lw  - f\Z2,  w € S^'m(A) 


It  has  been  shown  in  [10],  [4]  that  the  L^-projector  is  bounded  in  the 
uniform  norm,  provided  that  the  mesh  A is  quasiuniform.  The  boundedness 
of  the  discrete  L^-projector,  P^,  follows  as  in  Lemma  3.1.  Thus  we 
can  obtain  the  following  theorem,  referring  to  the  analysis  in  [ 5j  (Thm.  3.1 
there;  cf . also  [ 21 ] ) . 

Theorem  6.2. 

,(m+k+n)r 


Let  u € C' 


[ a,  b]  be  an  isolated  solution  for  (6. 1)-(1.  2) 


and  suppose  that  F(t,  zQ,  • • • > zm  ^ is  sufficiently  smooth  near  u. 

* ° q 

Then,  for  h small  enough,  the  approximate  solution  v e S (A),  solving 


(6.11) 


k+m' 

brn  , . m-1  . i2  i ^m  *,  , * ^m-1  * 

v(  • ) - F(  • , v,  . . . , D v)  I = | D v ( • ) - F(  • , v , . . . , D v 


v«Sr*  (A) 
k+m 


exists  and  is  unique  in  some  neighborhood  of  u,  and  Newton's  method 
(6.  5)-(6.6),  for  approximately  solving  (6.11),  converges  quadratically  in 


m 


some  neighborhood  of  v . We  have  D v = Y,  where  Y is  defined 

* 

by  (6.9).  Moreover,  the  error  estimates  (1 . 18)— (1 . 19)  hold  for  v here. 


7 . Numerical  examples 


The  collocation  method  of  de  Boor  - Swartz  [ S]  has  become  popular 

enough  in  recent  years,  so  that  we  feel  it  appropriate  to  use  this  method 

as  a frame  of  reference  in  comparing  the  discrete  least  squares  schemes. 

All  the  examples  here  are  of  second  order,  m = 2,  on  I = [0,1], 

with  Dirichlet  boundary  conditions 

(7.1)  u(0)  = u(  1 ) = 0 

and  are  smooth  enough  so  that  n = k. 

Unless  otherwise  specified  (as  is  in  example  2),  we  have  approximated 

using  quintics,  i.e.  k = 4,  and  on  a uniform  mesh.  Thus  h = 1/N, 

where  N is  the  number  of  intervals.  We  have  used  a B-spline  basis  [ 3] 

° ° q 

for  the  approximation  space  S = S.  (A): 

k+m 

5 = span{vl *’dim)  ' 

We  note  that  all  the  schemes  here  evaluate  the  differential  equation 

at  exactly  the  same  points,  namely,  the  Gaussian  points  { r . , } of  (1.9). 

The  collocation  solution  [ 5]  is  in  C*[0,1];  thus  it  has  enough  free 

parameters  to  satisfy  the  differential  equation  at  all  {t  },  i = 1 N, 

ij 

j = 1,  . . . , k.  The  resulting  matrix  of  the  linear  system  in  this  case  has 
an  "almost  block  diagonal"  structure,  elaborated  in  [7],  [6].  The 
method  proposed  in  [ 7 ] was  used  here  for  the  solution  of  the  linear  system 
of  algebraic  equations  resulting  from  the  collocation  scheme  C - 1. 

When  raising  the  continuity  at  the  knots,  by  requiring  the  approximate 
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* 


i 


1 


2 3 4 

solution  to  be  in  C [0,1],  C [0,1]  or  C [0,1],  we  obtain,  on  trying 
to  collocate  on  {t  },  an  overdetermined  system  of  linear  equations. 
This  system  is  solved  in  the  weighted  least  squares  sense,  as  explained 
in  the  introduction.  We  denote  the  schemes  obtained  by  C-2,  C-3  and 
C-4,  respectively.  The  normal  equations  were  solved  using  banded 
Cholesky  decomposition.  We  note  that,  compared  to  C-l,  the  condition 
number  of  the  matrices  in  the  other  schemes  has  been  essentially  squared, 
but  in  practice  we  have  not  sensed  any  effect  of  that.  This  indicates  a 
stable  structure  obtained  when  applying  B-splines  at  Gaussian  points. 

We  use  a notation  similar  to  [5]:  we  measure  the  quantities 

e(j)  = II DJ ( u - v*)||  ( ),  j = 0,1,2 

oc'  ' 

e^  = max  |dJ(u  - v*)(t)  I,  j = 0, 1,  2 
t«  A 


and  use  the  shorthand 


.5-4  for  .5  X 10 

(j) 


-4 


The  global  uniform  errors  e have  been  measured  on  an  "evaluation  grid' 
consisting  of  the  knots  (±  10  *°)  and  39  equally  spaced  points  between 
each  two  knots.  Further,  in  order  to  measure  the  rates  of  convergence, 
we  have  listed  in  row  i,  i > 1,  the  numbers 


V(j)  = log(e(ij)/e[J_)1)/log(N._1/N.) 

and  similarly  defined,  where  ef^  is  the  error  e^  and  N.  is 

A l l 
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the  number  of  (uniform)  intervals  N in  the  i-th  entry.  When  listing 
the  results,  we  have  omitted  the  column  e^  since  it  almost  always 
coincided  with  e^  i.e.,  the  largest  errors  in  the  second  derivatives 
were  usually  obtained  at  the  knots. 

Example  1. 

We  first  compute  for  the  simple  example  appearing  in  [17],  [5] 

(?•  2)  u"  - 4u  = 4 cosh(l) 

whose  solution  is 

(7-  3)  u(t)  = cosh(2t  - 1)  - cosh(l)  . 

Results  for  the  schemes  C-l  [5],  C-2,  C-3  and  C-4  are  contained  in 
Table  1,  where  under  "N"  we  have  listed  the  number  of  uniform  sub- 
intervals, and  under  "dim"  the  dimension  of  the  space  S?  ( 1)  (which 

k+m 

is  also  the  order  of  the  system  of  equations  solved).  The  actually 
computed  values  are 

aV  ’ • ‘ ’ “dim 

for  the  representation  of  the  approximate  solution 

...  dim 

(7-  4)  v'(t)  = Yj  «.?.( *)  • 

i = l 1 

The  theoretical  rates  of  convergence  and  superconvergence  (1. 18)-(1. 19) 
(and  also  (1. 11) - (1 . 12))  are  clearly  demonstrated  in  Table  1.  Note  also 
the  closeness  of  the  C-2  solution  to  the  collocation  solution  C-l.  They 
remained  close  in  all  the  examples  we  tried. 
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An  operation  count  for  the  collocation  and  discrete  least  squares 
methods  has  been  given  in  [ 18  ] . In  [6]  it  has  been  noted  that  the  collocation 
method  C-l  would  always  be  faster,  if  the  "almost  block  diagonal"  structure 
of  the  collocation  matrix  would  be  used  to  advantage  (rather  than  merely 
treating  the  matrix  as  "banded").  Our  timing  results  support  the  claims 
in  [ 18  ] , [ 6) : In  all  our  computations,  the  time  required  to  set  up  and 
solve  for  a certain  number  N moderately  increased  according  to  the 
following  order  from  left  to  right 

C-l,  C-4,  C-3,  C-2  . 

The  more  significant  measure  is,  however,  the  amount  of  work  required 

to  obtain  a required  accuracy.  This  measure  of  computing  time  would 

order  the  schemes  (according  to  all  our  computations)  in  the  following  sequence 

C-l,  C-2,  C-3,  C-4 

C-l  being  the  most  efficient,  and  C-4  the  least  efficient. 

The  above  observations,  plus  the  additional  flexibility  and  better  condi- 
tion number,  give  collocation  superiority  over  discrete  least  squares.  In  fact, 
the  only  reason  we  see  for  not  always  preferring  collocation  is  that  the  storage 
of  the  already  computed  solution  gets  more  expensive  as  the  dimension 
of  the  approximation  space  gets  higher.  The  C-2  solution  is  quite  close  to 
the  collocation  solution,  with  less  free  parameters,  and  still  partly  possessing 
the  superconvergence  phenomenon.  (Note,  though,  that  with  collocation  we 
can  get  superconvergence  also  for  the  m-th  derivative,  using  the  differential 
equation  and  the  values  of  the  approximate  solution  and  its  first  m - 1 


derivatives) . 


'v  . 


Table  1 


Errors,  convergence  rates  and  computing  time  for  example  1 . 


4 10 

5 12 

6 14 

7 16 

8 18 
9 20 
.0  22 
11  24 
.2  26 


Y(0) 

e<‘> 

Y(U 

e(2> 

. 35-5 

. 56-4 

. 29-2 

. 33-6 

5.8 

.81-5 

4.  8 

.61-3 

. 63-7 

5.8 

.20-5 

4.8 

. 20-3 

. 17-7 

5.9 

.68-6 

4.9 

. 84-4 

. 58-8 

5.9 

. 28-6 

4.9 

. 41-4 

. 23-8 

5.9 

. 13-6 

4.9 

. 23-4 

.11-8 

5.9 

. 68-7 

4.9 

. 1 3-4 

. 5 3-9 

5.9 

. 38-7 

4.9 

.84-5 

. 28-9 

5.9 

. 23-7 

4.9 

. 55-5 

. 35-5 

. 56-4 

.29-2 

. 34-6 

5.7 

.83-5 

4.7 

. 62-3 

. 64-7 

5.8 

. 21-5 

4.  8 

. 20-3 

. 17-7 

5.9 

.69-6 

4.9 

.85-4 

. 60-8 

5.9 

.28-6 

4.9 

. 42-4 

. 24-8 

5.9 

.13-6 

4.9 

.23-4 

. 11-8 

5.9 

. 69-7 

4.9 

. 14-4 

. 54-9 

5.9 

. 39-7 

4.9 

.85-5 

■ 29-9 

5.9 

. 23-7 

4.9 

. 56-5 

. 11-4 

.16-3 

. 55-2 

. 18-5 

4.  4 

. 35-4 

3.7 

.14-2 

. 43-6 

5.  1 

.11-4 

4.  2 

. 47-3 

.12-6 

5.  6 

. 39-5 

4.  5 

. 20-3 

. 45-7 

5.6 

. 16-5 

4.  7 

.99-4 

. 18-7 

5.8 

.78-6 

4.8 

. 55-4 

.8  3-8 

5.9 

. 41-6 

4.8 

. 33-4 

. 42-8 

5.9 

. 23-6 

4.8 

.21-4 

. 22-8 

5.8 

. 14-6 

4.  8 

. 14-4 

. 1 3-8 

5.8 

.89-7 

4.  8 

.94-5 

.78-9 

5.8 

. 58-7 

4.8 

. 67-5 

31-8 

11- 9  8.3 

12- 10  7.7 

19-11  8.2 

46-12  7.8 

13- 12  8.1 

46-13  7.9 

18-13  8.1 

79-14  7.8 

31-8 

77-10  9.1 

81-11  7.8 

15-11  7.6 

39-12  7.5 

11-12  7.8 

41-13  7.7 

18-13  7.1 

61-14  10.2 


78-7 

33-8 

33-9 

57-10 

13-10 

39-11 

13-11 

52-12 

23-12 


“ 6 . 4 | 
5.0 
5.  3 | 
5.  5 
5.  6 
5.6 
5.  7 
5.  7 


a* 


Table  1 (cont. ) 


N dim 


2 

5 

3 

6 

4 

7 

5 

8 

6 

9 

7 

10 

8 

1 1 

9 

12 

10 

13 

1 1 

14 

12 

1 5 

13 

16 

14 

17 

15 

18 

e<°> 

V(0) 

e<» 

(1) 

\ 

(2)  (2) 
e v 

si01  vi0) 

~ n 

(i)  u; 

6 A VA 

.11-4 
.47-5  2.1 

.71-6  6.6 

.26-6  4.6 

. 16-3 

.69-4  2.0 

.16-4  5.1 

.60-5  4.4 

. 55-2 
.11-2  3.9 

.41-3  3.5 

.16-3  4.2 

. 11-4 

.31-5  3.1 

.71-6  5.2 

.17-6  6.5 

j 

. 58-6 

.21-4  -8.8 

.41-5  5.6 

.20-5  3.2 

. 80-7 

6.4 

. 24-5 

5.0 

.81-4  3.8 

.60-7  5.6 

.74-6  5.6 

. 36-7 
.16-7 

5.  2 

6.  1 

.12-5 
. 62-6 

4.7 

4.8 

.46-4  3.8 

.27-4  3.9 

.23-7  6.3 

.11-7  5.5 

.38-6  4.3 

.20-6  4.8 

00 
ro 

1 

oo 

5.  6 

.35-6  4.8 

.18-4  3.8 

.53-8  6.0 

.12-6  4.6 

. 44-8 

5.9 

.21-6  4.8 

.12-4  3.8 

.30-8  5.6 

.71-7  4.8 

. 25-8 

5.8 

.13-6  4.8 

.81-5  3.8 

.17-8  5.9 

.45-7  4.7 

. 1 5-8 

5.9 

. 88-7 

4.8 

.58-5  3.9 

.10-8  5.7 

.30-7  4.8 

.9  5-9 

5.8 

. 60-7 

4.8 

.43-5  3.9 

.64-9  5.8 

.20-7  4.8 

. 61-9 

5.9 

. 42-7 

4.9 

.32-5  3.9 

.42-9  5.8 

.14-7  4.8 

. 41-9 

5.9 

. 30-7 

4.9 

.25-5  3.9 

.28-9  5.8 

.10-7  4.8 

Example  2 . 


Here  we  consider  a less  trivial  problem  of  boundary  layer,  similar 
to  that  considered  in  [ 17  ] , [ 20  ] 

(7.5)  eu"  - (2  - t2)u  = -1  + j(Z  ~ t2)(l  - t)  + (1  - t2)exp{-  ^ } 

with  a very  good  approximation  to  the  exact  solution,  for  a small  e > 0, 
given  by 

(7.6)  u(t)  = — l—  - exp{-  + j(t  - 1)  . 

The  solution  varies  very  fast  between  t = 1 - 0(*/e)  and  t = 1.  We 
- 8 

took  e = 10  . (Note  that  the  "trouble  making"  term  of  u'  is  absent 

in  (7.  5).  However,  it  is  still  of  interest  to  observe  the  performance  of 
the  schemes  for  this  rapidly-varying  solution). 

We  followed  [17]  by  choosing  the  highly  nonuniform  partition  (N  = 23) 
0,  .1,  .3,  .5,  .7,  .8,  .9,  .95,  .97,  .98,  .985,  .99,  .9925,  .995, 

.996,  .997,  .9975,  .998,  .9985,  .999,  .9995,  .9997,  .9999,  1.0 
and  have  computed  piecewise  cubic  (k  = 2)  and  quintic  (k  = 4) 
approximate  solutions  by  the  collocation  and  discrete  least  squares  schemes 
on  Gaussian  points.  Results  are  listed  in  Table  2.  Note  that  we  may  already 
be  in  the  asymptotic  range  of  superconvergence  for  "quintic  C-l" 
in  e ^ and  and  for  "quintic  C-2"  in  Note  also  that 

the  errors  listed  in  Table  2 are  absolute,  rather  than  relative,  and  that  we 
have  for  t » 1 - *^e 


u = 0(1),  u'  = 0(1  /S/e),  u"  = 0(l/e)  . 


Table  2 


Boundary  layer  problem,  example  2 (N  = 23)  (e  = 10  ) 


method 

dim 

.<*> 

e<2> 

e<3> 

e(°) 

eA 

4U 

cubic  C-l 

46 

. 39-2 

.95+2 

. 67+7 

. 19-2 

. 27+2 

cubic  C-2 

24 

. 69-2 

. 14  + 3 

. 67+7 

. 31-2 

.8  5+2 

quintic  C-l 

92 

. 29-4 

. 12  + 1 

. 17+6 

. 63-5 

. 62-1 

quintic  C-2 

70 

. 36-4 

. 14  + 1 

. 1 3+6 

.93-5 

.72+0 

quintic  C-  3 

48 

.15-3 

. 39+1 

. 21+6 

.11-3 

. 32  + 1 

quintic  C-4 

26 

. 94-3 

. 18+2 

.72+6 

.92-3 

.93  + 1 

Example  3. 

We  now  turn  to  a physically  interesting  nonlinear  problem  [17],  [14],  [9] 

(7.7)  u"  + t exp{u}  = 0,  t > 0 . 

The  solutions  of  this  problem  are  given  by 

(7. 8)  u(t)  = -2  log{cosh[  (x  ~ ]/cosh[~^  ] } 

where 

(7.9)  0 = iv/2t  cosh(-^)  . 

There  exists  a critical  value  of  t,  such  that  equation  (7.9)  has 

two,  one  or  no  solutions  according  to  t being  t < t , t = t or 

c c 

t > t , respectively.  The  above  solutions  are  all  the  possible  solutions 

of  (7. 7)-(7. 1)  [ 9] ; i.  e. , beyond  there  is  no  solution  to  the  problem. 

For  the  two  solutions  u,,  u which  do  exist  when  t < t , we  have 

1 z c 

(7.10)  u](t)<u2(t),  o < t < 1 

and  u^,  is  unstable  [13]  . The  numerical  solution  will  converge  to  u^. 
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Clearly,  the  more  we  increase  t,  t < t 


the  more  difficult  it 


becomes  to  obtain  an  accurate  numerical  solution.  In  order  to  find  the 
critical  value  t^,  let  0c  be  the  corresponding  solution  of  (7.9). 
We  have 

“ 0 = ctgh(^  0 ),  (ctgh  = ?OS/1 ) 

4c  4c’  sinh 


t = 8/sinh2(4  0 ) . 
c 4 c 


The  resulting  values  are,  approximately, 


t a 3.  51383,  0 * 4.8. 

c c 

We  have  computed  with  two  values  of  t ; 

(i)  t = 1 (then  for  u,  0 = 1.5171645990508) 

(ii)  t = 3.  5 (then  for  u , 0 = 4.5518536628383)  . 

We  have  applied  Newton's  method,  starting  with  vQ  £ 0 and 
computing  v j the  approximate  solution,  according  to  the  appropriate 

scheme  C - i,  of  the  linear  problem 


(7.10) 


y"  + (t  exp{vr))y  = (vf  - 1)t  exp{vf} 
y(0)  = y(l)  = 0 . 


The  iteration  process  was  stopped  when  v agreed  with  v to  within 

r+i  r 

at  least  13  significant  decimal  digits. 

(i)  For  t = 1,  results  are  accumulated  in  Table  3.  In  all  cases  it 
took  3 iterations  for  the  sequence  {v^}  to  converge.  The  properties  of 
convergence  and  superconvergence  are  demonstrated. 


Table  2 


Boundary  layer  problem,  example  2 


(N  = 23)  (e  = 10~8) 


method 

dim 

e“> 

e(2) 

e(3» 

40) 

4l) 

cubic  C-l 

46 

. 39-2 

.9  5+2 

. 67+7 

. 19-2 

. 27+2 

cubic  C-2 

24 

. 69-2 

. 14  + 3 

. 67+7 

. 31-2 

.85+2 

quintic  C-l 

92 

. 29-4 

. 12+1 

.17+6 

. 63-5 

. 62-1 

quintic  C-2 

70 

. 36-4 

. 14  + 1 

. 1 3+6 

.93-5 

.72+0 

quintic  C-3 

48 

.15-3 

. 39+1 

. 21+6 

.11-3 

. 32+1 

quintic  C-4 

26 

.94-3 

. 18+2 

.72+6 

.92-3 

.93  + 1 

Example  3. 

We  now  turn  to  a physically  interesting  nonlinear  problem  [17],  [14],  [9] 

(7.7)  u"  + t exp{u}  = 0,  t > 0 . 

The  solutions  of  this  problem  are  given  by 

(7.8)  u(t)  = -2  log {cosh[  (x  - j )|  ]/cosh[  ~ ] } 
where 

(7.9)  9 = 'v/2t  cosh(-^)  . 

There  exists  a critical  value  of  t,  t such  that  equation  (7.9)  has 

two,  one  or  no  solutions  according  to  t being  t < t t = or 

t > t , respectively.  The  above  solutions  are  all  the  possible  solutions 

of  (7 . 7 )— (7 . 1)  [ 9] ; i.  e. , beyond  there  is  no  solution  to  the  problem. 

For  the  two  solutions  u,,  u_  which  do  exist  when  t < t we  have 

1’  2 c’ 

(7.10)  Uj(t)  < u2(t),  0 < t < 1 

and  u^  is  unstable  [13].  The  numerical  solution  will  converge  to  u^. 


-41- 


(ii)  For  r = 3.  5,  results  are  listed  in  Table  4.  As  expected,  the 
problem  is  harder  here;  nevertheless,  the  schemes  still  behave  well.  The 
number  of  iterations  needed  for  convergence  of  the  Newton  process  is 
considerably  higher  here,  and  is  listed  under  the  heading  "iter". 

The  above  computations  were  performed  in  double  precision  on  ihe 
Univac  1110  computer  at  the  University  of  Wisconsin-Madison. 

The  author  wishes  to  thank  Professor  C.  de  Boor  for  many  helpful 
discussions. 
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Table  3 


Nonlinear  example,  t = 1. 


N dim 


4 1 6 
C-l  6 24 
8 32 

4 13 
C-2  6 19 

8 25 

4 10 
C- 3 6 14 

8 18 


43-8 

39- 9  5.9 

71- 10  6.0 

46-8 

40- 9  6.0 

72- 10  6.0 


(1)  (1) 

e \ 


1 4-6 

19-7  4.9 

45-8  5.0 

14-6 

19-7  5.0 

45-8  5.0 

83-6 

13-6  4.5 

34-7  4.8 

12-7  4.9 


(2)  (2) 

e \ 


1 3-4 

27-5  3.9 

87-6  4.0 


65-11, 

23-12  8.2 

16-13  9.4 


48-11 

22-12  7.6 

48-1  3 5.3 


. 59-11 

. 44-7 

4.0 

. 22-12 

8.  1 

. 43-8 

5.7 

4.  0 

.15-13 

9.  3 

.77-9 

5.9 

. 33-7 

. 29-6 

3.  6 

. 36-8 

5.  4 

. 37-7 

5.  1 

3.  8 

. 69-9 

5.7 

.76-8 

5.  5 

3.9 

. 19-9 

5.  8 

. 22-8 

5.7 

.7  3-7 

.75-6 

4.4 

. 47-8 

6.8 

.71-7 

5.  8 

4.0 

. 85-9 

6.  0 

.13-7 

5.  8 

4.  1 

. 21-9 

6.  3 

. 36-8 

5.8 

4.0 

.70-10 

6.0 

.12-8 

5.  8 

4.  1 

. 27-10 

6.  1 

. 51-9 

5.  8 
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